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Abstract 


We study the thermal boundary conduction in one-dimensional harmonic and b 4 lattices, both 
of which consist of two segments coupled by a harmonic interaction. For the ballistic interfacial 
heat transport through the harmonic lattice, we use both theoretical calculation and molecular 
dynamics simulation to study the heat flux and temperature jump at the interface as to gain 
insights of the Kapitza resistance at the atomic scale. In the weak coupling regime, the heat 
current is proportional to the square of the coupling strength for the harmonic model as well as 
anharmonic models. Interestingly, there exists a negative temperature jump between the interfacial 
particles in particular parameter regimes. A nonlinear response of the boundary temperature jump 
to the externally applied temperature difference in the b 4 lattice is observed. To understand the 
anomalous result, we then extend our studies to a model in which the interface is represented by 
a relatively small segment with gradually changing spring constants, and find that the negative 
temperature jump still exist. Finally, we show that the local velocity distribution at the interface 
is so close to the Gaussian distribution that the existence/absence of local equilibrium state seems 
unable to determine by numerics in this way. 
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I. INTRODUCTION 


fl, 


Since the pioneering observation between liquid helium and a metal [lj], thermal bound¬ 
ary resistance, namely, Kapitza resistance, has been extensively studied theoretically and 


experimentally Q . With the rapid development of modern electronic technology, there has 
been much need and interest in understanding the fundamental nature of thermal boundary 
conductance since it has been a significant obstacle in designing the micro- or nano- scale 
electronic chips. 


mismatch model 


wo phenomenological models, the acoustic mismatch model [3| and diffuse 


a, 


have been proposed to study the mechanism of the thermal boundary 
conductance. However, due to the neglecting of atomic details at the interface, they both 
offer limited accuracy, particularly, for nanoscale interfacial resistance [i]. To understand 


the mechanism of thermal boundary conductance at the atomic 


evel, many studies have 
6]. Most of the previous 


been done in one-dimensional lattices via different methods m 
studies focus on the effect of the interface on the steady-state heat flux and little attention 
has been paid to the temperature jump between the interface from the atomic viewpoint. 


On the other hand, heat conduction in low-dimensional dynamical systems has become 


L 
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re subject of a large number of theoretical and experimental studies in recent years 


EL 


9]. An exact approach to interacting Hamiltonian systems is so far unavailable except 
for harmonic crystals. A meaningful definition of local temperature depends on the local 
thermal equilibrium and it is difficult to give a microscopic derivation of the condition 
in general jlO|. With the usual definition of local temperature i.e., the mean local kinetic 
energy, the temperature profile may show some unexpected features, such as the temperature 
oscillations in the steady state of alternate mass hard particle gas jll|. in the Fermi-Pasta- 


Ulam chain 


12] and in harmonic chain with alternating mass [13 ]. 


In the present study, we study the heat flux and temperature jump at the interface as to 
gain insights of the Kapitza resistance at the atomic scale via theoretical calculations and 
molecular dynamics simulations. We fold that there exists a negative temperature jump 
between the interfacial particles in particular parameter regimes. A nonlinear response of 
the boundary temperature jump to the externally applied temperature difference in the 0 4 
lattice is observed. Note that, although the interface between two segments is not well 
defined in one-dimensional Hamiltonian systems, our studies can give some insights into the 
thermal boundary resistance in real systems. 
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The paper is organized as follows. In Sec. [T]] we define the model and give the methods 
for theoretical calculations and molecular dynamics simulations. In Sec. Illll we demonstrate 
the existence of negative temperature jump in both the harmonic and 0 4 model. Finally, we 
give a brief summary and discussion in Sec. IIV1 


II. MODEL AND METHODS 


We study the non-equilibrium steady state of a one dimensional chain consisting of two 
coupled lattice, 


2 

H = H l + H r + —k c (xjy/2 — %N/2+l) ■ 


( 1 ) 


The Hamiltonian for the left and right segments are given by 

JV/2-1 


i =1 v 7 i=0 


Y (Xi-X i+1 f, 


( 2 ) 


N 


Hr = 


l Pi , f R 2 , 


2 m 


N 


[Xi 


Xi+iY 


i=N/ 2+1 x 7 i=N/ 2+1 

where Xi denotes the displacement of the i-th particle from its equilibrium position. Fixed 
boundary conditions are taken, i.e., x 0 = xjv + i = 0. The particle 1 and N at the two ends 
are connected to the heat baths at temperature T L and T R , respectively. The heat baths 
are modeled by the Langevin equations corresponding to Ohmic baths, i.e., the self energy 
of the baths are E(cu) = iya; 8]. 

When Xl = X R = 0, the on-site potential and inter-particle interaction are all quadratic. 
In the classical limitation, the steady heat current from left to right reservoir can be given 
by the Langevin equations and Green’s function (LEGF) method si. [ldj] 


J = 


k B (T L -T R ) 


7T 


duTi 


with 


\_—u 2 Ms + K s - - E±((n)] ’ 

r L,fl(w) = Im (E+ R (a;)) , 


( 3 ) 

( 4 ) 

( 5 ) 


where Ms and Kg denote the mass matrix and force constant matrix of the system. Note 


that G s , E l r are all N x N matrices for one-dimensional systems. The only non-zero 
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element of E± R are respectively [Ef]^ = E = i^/ui and [E^]jv,at = E = i'yu. 7 is the 
coupling strength of the first and N -th particle to the left and right reservoirs, respectively. 
The velocity-velocity correlation and position-velocity correlation are: 


K = (X S X J) = 


^bTl 


7T 


duuGi; (o;)r l(uj)G s (a;) 


k bTr 


7T 


duuGg (a;)T r(u)G s ( u ), 


( 6 ) 


C = <■ X S X J) = 


ik B T L 


7T 


duGg (a;)T l{u)G s ( u ) 


( 7 ) 


+ 


ik R T R 


7T 


duGg(cu)T r(oj)G s ( u ). 


The correlation function K can be used to define local energy density which can in turn be 
used to define the local temperature, i.e., 


Tj = mK, 


l,l J 


( 8 ) 


and C gives the local heat current density 


15 


18]. We integrate Eq. fl3]) and Eq. (J 6 J) numeri¬ 


cally to obtain the steady state heat current and local temperature, for which the rectangular 
method is used 19]. We also verify that the local current, obtained by integrating Eq. (IT|) 
numerically, is the same along the chain, which is one of the properties in the steady state. 

When Al, Xr 7 ^ 0, we apply the non-equilibrium molecular dynamics simulation(NEMD) 
to the system, for which the Langevin heat baths are used at the two ends of the chain. The 
equations of motion are given by 

dH 


mxi = 


dxi 


7 i x i + du 


( 9 ) 


where 7 j = 7 (^ 1 ,* + lv,i) and r/, ; = tjlSi ^ + r] R 5N^. The noise terms r] L R denotes a Gaussian 
white noise with zero mean and variance of 27 k R Ti R . The local heat flux is given by 
j = (F(x i+ i — Xi)v i+ 1 ), where F(x) = —V'{x) and the notion (...) denotes a steady-state 
average. The equations of motion (Eq. (j9])) are integrated by using a second-order Stochastic 
Runge-Kutta algorithm 20]. At steady states, the numerically computed local heat flux is 
always constant along the chain, and the local temperature is defined as T* = m(x |). To 
compute the boundary temperature jump, i.e., A T b = T N / 2 — T/v/ 2 + 1 , the relaxation and 
average time must be both long enough. In what follows we set m — 1, k R — 1, and 7 = 1 . 
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FIG. 1. Heat flux as a function of interfacial coupling k c via the LEGF approach as k c approaches 
to zero. Symmetry: /cl = kR = 1,/l = f R = 0; Asymmetry: /cl = l,k R = 2 ,f L = f R = 0; 
Symmetry with on-site potential: /cl = 1, k R = 1, fL = f R = 2. Asymmetry with on-site potential: 
/cl = 1, k R = 2, /l = f R = 2. For all cases, we set Tl = 2,T R = 1 and N = 64. 

III. RESULTS AND DISCUSSIONS 


For many devices of several segments, interfacial coupling is pretty weak, which indicates 
that k c is far less than ki and k r in our model. So it is desirable to study the thermal 


transport through atomic chains in the weak coupling regime. It has been shown that 21] 


the heat current is proportional to the square of the coupling strength in one-dimensional 
weakly-coupled chain with the Morse on-site potential by a phenomenological analysis. Does 
this square-law relation between heat current and coupling strength is still valid in the weak 
coupling limit when the anharmonic on-site potential is absent? As shown in FigJU we 
plot the heat current as a function of the coupling strength in the weak coupling limit by 
integrating Eq. (|3]) numerically. It turns out that the square law relation is still valid when 
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the anharmonicity is absent. Further more, the square law relation still holds when the 
system consists of symmetrical/asymmetrical segments with/without an on-site potential. 
Fig. [2] shows the steady-state heat current and the boundary temperature jump as a function 



c c 


FIG. 2. The heat flux (left) and temperature jump between the IV/2-th and (N /2 + l)-th particles 
(right) as a function of k c via both the LEGF approach and MD simulations. Here /l = Ir = 
0, Tl = 2 ,Tr = 1, Xl = \r = 0, and N = 64. 

of the coupling strength k c . The reason to carry out both theoretical calculation and NEMD 
simulation is to verify that the results we obtained are from physical reasons rather than 
uncertain numerical reasons because the temperature jump between the N/ 2-th particle 
and the (N /2 + l)-th particle requires a highly accurately performed simulation for its 
sensitive to heat fluctuation when k c approaches to kR. By inspecting the figure, we can 
see that theoretical calculations and MD simulations agree well with each other. The heat 
current increases at first, then arrives at a maximum value, and then slightly decreases 
with the increase of k c . As depicted in Fig. [2] the maximum heat current occurs at k c = 
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2 kLkji/(kL + kji), which agrees with the result obtained in [5j by the scattering boundary 
method. Furthermore, both theoretical calculations and MD simulations indicate that there 
is a negative boundary temperature jump, i.e., A T b < 0, when k c approaches to k r. The 
word “negative” is in contrast with “normal” heat conduction that the direction of the heat 
flow is from hot to cold regions. In fact, similar negative temperature jumps occurs in 
several systems, for example, the temperature jump between the second and third particle, 
and between (TV — 2)-th and (TV — l)-th particle in the uniform harmonic chain jj] coupled 
with reservoirs, and the temperature oscillations in the steady state of hard particle gas 
111 ], the Fermi-Pasta-Ulam chain 12], the harmonic chain 13] with alternating mass. To 



FIG. 3. The contribution of normal modes to local temperature at the interface. Here ki, = kn = 1, 
fh = fn = 0, Xl = Xr = 0, Tl = 2, Tr = 1 and TV = 64. 

understand the negative temperature jump at the interface, we need to inspect the concept 
of the local temperature further. The local temperature of the i-th particle can be written 
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as Ti = Ai(uj max ), with 

AiH = 2 doj' (^^u'G%{u/)T l {u/)Gs{u/) + ^^uV%{u')T r (uj')G- s {uj')^ , (10) 

and Umax is the top boundary of the phonon spectra. The kinetic energy of a particle gets 
contributions from all the modes, and the net result depends on the distribution of energy 
in the different modes. As shown in Fig. [3] we plot the contribution of normal modes to 
the local temperature for the (AT/2)-th and {N / 2 + l)-th particles by integrating Eq. (TTOl) 
numerically. As we can see, “equipartition” among phonon modes, i.e., each normal mode 
shares the same average kinetic energy, is not satisfied for k c = 0.5 and k c = 1 shown by the 
nonlinear behaviors of Aj(cu)(i = N/2, N/2 + 1) in the high frequencies region. Surprisingly, 
for the case of k c = 1.5, A*(a;) exhibit almost linear behavior with the increasing of ui, 
indicating that the contribution to the local temperature from possible phonon modes are 
closely equivalent. Comparing with k c = 0.5, we can see that the high frequency normal 
modes are suppressed more dramatically than the low frequency normal modes for k c = 1.5 
and k c = 2, and the turning of Aat/ 2 and Ajv/ 2 +i in the high frequencies region indicates the 
negative temperature jump between the AT/2-th and {N/2 + l)-th particle. 

It would be interesting to see if the negative temperature jump is an artificial effect due 
to the integrability of the harmonic system. Thus we conduct similar studies in the </> 4 
lattice, which has additional nonlinear on-site potential on each site in comparison with 
the harmonic system. We plot the boundary temperature jump AX), as a function of the 
external temperature difference AT = Tl — T R and some typical temperature profiles in 
Fig. SI One can see that the boundary temperature jump is proportional to AT for k c < 1 
and proportional to {—AT) for k c > 1 when AT is small, which are typical linear-response 
behaviors shown in harmonic models. With the increasing of AT, the linear behavior of 
AX), no longer holds for the </> 4 lattice. Note that negative temperature jump occurs when 
k c > 1.3 and the absolute value of AT h nonlinearly increases as AT increases. 

So far our discussions is based on the model consisting of two segments with a single 
harmonic coupling, which inevitably leads to the argument that the origin of negative tem¬ 
perature jump comes from the ill-defined interface of the two-segment model with a sharp 
discontinuity of the interfacial coupling. In what follows we propose an extended model to 
show that it is not the case. Actually, in a practical consideration, the interface may be a 
junction which is small compared with the two sub-lattices. So we divide our system into 
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FIG. 4. ATj, as a function of AT for different k c in the iji 4 lattice , with Tl = Tr + AT, Tr = 
0.5, Al = = 1, /cl = = 1 and N = 64. 

three regions, say, two sub-lattices and a junction. The particle number of the junction is 
small compared with the two sub-lattices. The spring constant of the intermediate segment 
varies smoothly, which is done by setting the spring constants of the intermediate junction 
by ki = exp (— (i — N/ 2) 2 /50) + 1, where i represents the index of particles. The NEMD 
simulation results of harmonic and 0 4 lattices are presented in Fig. [5] and Fig. [6j respectively. 
As we can see, negative temperature gradient still exists for both harmonic and </> 4 lattices 
within the interfacial segment. 

As mentioned above, a meaningful local temperature can be defined only in systems 
exhibiting local thermal equilibrium. And we known that, if the system can exhibit local 
thermal equilibrium, the local distribution should be gaussian and all even moments can be 
obtained based on the second moment. We can then use T- 2 ' 1 = m{xf), = m ((x; 4 )/3) 1//2 , 
and t/ 6) = m ((i; 6 )/15) 1 ^ 3 to define local temperature equivalently. So we carry out NEMD 
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FIG. 5. The temperature profile of the harmonic chain with an interfacial junction, whose spring 
constants are smoothly varied. The distribution of spring constants for the whole system is given 
as follows: ki = ki = 1 for 1 < z < 7IV/16; ki = exp (— (i — N/ 2) 2 /50) + 1 for 71V/16 < i < 9IV/16; 
and ki = kn = 1 for 9IV/16 < i < (N — 1), where i is the index of particle number and N = 256. 
The first three even moments of velocity are given by = m{xf) for the second moment, 
= m ((xj 4 )/3)for the fourth moment and T = m ((xj 6 )/15)for the sixth moment, 
respectively. Here Tl = 2,Tr = 1, and N = 256. 

simulation for both the harmonic and </> 4 lattice, and plot the local temperature defined 
by the first three even moments of velocity, namely, the second, fourth and sixth moment 
in Fig. ED and Fig. |HJ To our surprise, the local temperature defined by T^ 2 \ Td) and 
at the boundary particles agree well with each other. The deviation at the interface 
is not significant in comparison with that inside segments. The result indicates the local 
distribution is or at least very close to the Gaussian, which cannot be well distinguished by 
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FIG. 6. The temperature profile of the ^> 4 lattice with an interfacial junction, whose spring constants 
are smoothly varied. The distribution of the spring constants in the interfacial junction is the same 
as that for Fig. [5j Here = Xr = 1 ,Tl = 2,Tr = 1, and N = 256. 

numerics and should recourse to more careful theoretical studies of local distribution in the 
future. 


IV. SUMMARY 

We have studied interfacial thermal conductance in one-dimensional inhomogeneous sys¬ 
tems by using both theoretical calculations and MD simulations. In the weak coupling limit, 
theoretical calculations show that the heat current is proportional to the square of the cou¬ 
pling strength in the absence of anharmonicity. A negative temperature jump between the 
interfacial particles occurs in both the harmonic and 0 4 lattices. As to understand the coun- 
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terintuitive observation, we have investigated the contribution of normal modes to the local 
temperature at the interface. It is shown that the high frequency modes make dominant 
contribution when the coupling strength is small, however, the contribution of each mode is 
almost equivalent when the coupling strength is strong. We have confirmed the occurring 
of the negative temperature jump is not trivially artificial due to the integrability of the 
system or the sharp discontinuity of the interfacial coupling by extending the system to a 
model consisting of two sub-lattices and an intermediate junction for both the harmonic and 
cj) A lattices. 


One should reexamine the notion of temperature as to understand the anomalous negative 
temperature jump, which seemingly indicats that heat flows against a local temperature gra¬ 
dient in a small scale. On the one hand, from the viewpoint of traditional thermodynamics, 
local temperature should be defined in a “cell” which should be macroscopically infinitesimal 
but contain enough microscopic degrees of freedom. Such kind of cell is, strictly speaking, 
not well defined for our microscopic model due to the large atomic-scale fluctuations and 
the word ’’local” defined for a single oscillator loses its inherent meaning. On the other 
hand, we stress that the traditional definition of local temperature with respect to the ki¬ 
netic energy of an oscillator is still in the framework of equilibrium thermodynamics. The 
anomalous phenomenon may partly comes from the definition as used here, which lacks a 
complete description of the nonequilibrium steady state. A new definition of “nonequilib¬ 


rium temperature” might be taken into consideration on this count 


22 


23], especially when 


one take notice of the temperature profile for the middle region of the intermediate junction, 
which is anomalously smaller than T ri as shown in Fig. [5] However, whether the concept of 
(local) temperature can be extrapolated beyond local equilibrium or should be modified in 
the nonequilibrium systems is still an open question. 
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